Universal mean moment rate profiles of earthquake ruptures 
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Earthquake phenomenology exhibits a number of power law distributions including the Gutenberg- 
Richter frequency-size statistics and the Omori law for aftershock decay rates. In search for a basic 
model that renders correct predictions on long spatio-temporal scales, we discuss results associated 
with a heterogeneous fault with long range stress-transfer interactions. To better understand earth- 
quake dynamics we focus on faults with Gutenberg-Richter like earthquake statistics and develop 
two universal scaling functions as a stronger test of the theory against observations than mere scal- 
ing exponents that have large error bars. Universal shape profiles contain crucial information on 
the underlying dynamics in a variety of systems. As in magnetic systems, we find that our analysis 
for earthquakes provides a good overall agreement between theory and observations, but with a 
potential discrepancy in one particular universal scaling function for moment-rates. The results 
point to the existence of deep connections between the physics of avalanches in different systems. 
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INTRODUCTION 

Earthquake phenomenology is characterized by several 
power law distributions. The most famous of these are 
the frequency-size distributions (i.e. histograms) of re- 
gional and global earthquakes 0, |||, and the modified 
Omori law for the aftershock decay rate around large 
rupture zones 0, Q . Using the seismic moment Mq for 
the earthquake size, the frequency-size distributions (or 
moment histogram) has the form (e.g., Q) 



n(M ) ~ Mi 



(1) 



where Mq ~ J^. AttjAA; with Auj and AAi being the 
local slip and rupture area during an earthquake, respec- 
tively. A related, more commonly used form in terms of 
earthquake magnitude M is 



log(n(M)) = a - bM 



(2) 



where n(Mo)dMo = n(M)dM, the constant a character- 
izes the overall rate of activity in a region, and the "&- 
value" gives the relative rates of events in different mag- 
nitude ranges. Using the observed moment-magnitude 
scaling relation M ~ 2/3log(Mo) for large earthquakes 
[3- 0- E3| - the exponent (3 of jtj is related to the b- value of 
@ as b = 1.5/3. The modified Omori law for aftershock 
decay rates is: 



AN/ At ~ K/{t + cf 



(3) 



where N is the cumulative number of aftershocks, t is the 
time after the mainshock, and K, c, and p are empirical 
constants. The exponents in (QJ and |J3J| are stable for 
data collected over large space-time domains, with some 



clear deviations from global averages related to fault- 
ing type and regional properties For example, the 
6-values of strike-slip, thrust, and normal earthquakes 
with depth < 50 km in the global Harvard catalogue are 
about 0.75, 0.85, and 1.05, respectively Q. As another 
example, regions with high heat flow often have short af- 
tershock sequences with relatively large exponent (e.g., 
p > 1.25), while regions with low heat flow have long 
aftershock sequences with low exponent (e.g., p < 0.9) 
0,0] ■ The association of earthquake statistics with power 
law relations like Eq. Q and Eq. © led some to sug- 
gest that earthquake dynamics is associated with an un- 
derlying critical point [1, El El 0] • However, power 
law distributions can be generated by many other mecha- 
nisms |l5l Il6j and it is important to develop criteria that 
can provide stronger evidence for or against the associa- 
tion of earthquakes with criticality. 

Recently, enough data have been collected to extract 
statistics of earthquakes on individual fault zones occupy- 
ing long (order 100 km) and narrow (order 10 km) regions 
of space. Wesnousky and collaborators 0, [w| found 
that the frequency-size statistics of earthquakes on highly 
irregular fault zones, with many offsets and branches, as 
the San Jacinto fault in California, are also described 
by the Gutenberg-Richter power law relation up to the 
largest events. However, relatively regular fault zones 
(presumably generated progressively with increasing ac- 
cumulated slip over time), such as the San Andreas fault 
in California, display power-law frequency-size statistics 
only for small events. These occur in the time intervals 
between roughly quasi-periodic earthquakes of a much 
larger "characteristic" size that is related to large-scale 
dimensions of the fault zone 

II 13 IH G3 (If the ratio 
of the mean divided by the standard deviation of the dis- 
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tribution of time intervals between characteristic earth- 
quakes is larger than 1, the distribution is referred to as 
quasi-periodic |Ti|). Earthquakes of intermediate mag- 
nitude are typically not observed on these faults (other 
than, perhaps, during aftershock sequences). The corre- 
sponding frequency size statistics are called the "charac- 
teristic earthquake" distribution 0, 0] . 

Previously these two types of behavior on individ- 
ual fault zones have been modeled as statistics close-to 
and far-from an underlying critical point 0, 0] , using a 
model for a strike-slip fault that incorporates long-range 
interactions and strong heterogeneities H3,E3| The dif- 
ferent dynamic regimes were associated with a compe- 
tition between failure-promoting effects of elastic stress- 
transfer or dynamic weakening, and the opposing effect 
of strength inhomogcneitics in the fault structure. Fisher 
et al. y found that near the critical point the frequency- 
size statistics follow a power law distribution (with a cut- 
off at large magnitude), with the same scaling exponent 
of observed data for strike-slip faults (i.e., a b- value of 
0.75). A similar form of frequency-size statistics and pre- 
dicted b- value were obtained also for a critical parameter 
value in a stochastic branching model [2f)j |. 

To provide an improved understanding of earth- 
quake dynamics that can suggest additional observables, 
we focus on faults with Gutenberg-Richter like earth- 
quake statistics (i.e. near-critical behavior) and develop 
two universal scaling functions associated with mean 
moment-rate time profiles at either fixed total moment or 
fixed total earthquake duration. Universal scaling func- 
tions (or shape profiles) give important information on 
the underlying dynamics, and may be found in solar 
flares in astrophysics |2l| . price fluctuations in financial 
markets j^, Barkhausen noise in magnets j2^|, and, as 
shown here, also in earthquake phenomenology. If the 
behavior of fault zones with earthquakes following the 
Gutenberg-Richter statistics is indeed critical, then the 
shapes of these functions should be as universal as the 
exponent [3 in Eq. JQ. Comparing the scaling functions 
in our earthquake model to observations constitutes a 
much stronger test of the theory than merely comparing 
a discrete, finite set of critical exponents. 

In this study we compute the scaling functions for both 
model predictions and observational data and compare 
the results. In section D of the paper we review the model. 
In sections D to D we extend the model, while also introduc- 
ing an extended phase diagram for the model dynamics, 
and the scaling behavior on long length scales. In section 
B we introduce the universal scaling functions and their 
scaling forms, and in section B we extract the functions 
from both simulation and observational data. Finally, in 
section D we discuss the results and emerging new ques- 
tions. 



EARTHQUAKE MODEL 

The model we use was developed originally by Ben- 
Zion and Rice 0,0], who suggested that a narrow irreg- 
ular strike-slip fault zone of horizontal length L and verti- 
cal depth W may be represented by an array of N ~ LW 
cells in a two dimensional plane, with constitutive pa- 
rameters that vary from cell to cell to model the disorder 
(offsets etc.) of the fault zone structure (FIG.nj. The 
cells represent brittle patches on the interface between 
two tectonic blocks that move with slow transverse ve- 
locity v in the x direction at a great distance from the 
fault. The interaction between cells during slip events is 
governed by 3-D elasticity and falls off with a distance 
r from the failure zone as -3-. These interactions are 
sufficiently long range that scaling in mean field theory 
(where the interaction range is set to infinity) becomes 
exact, up to logarithmic corrections, in the physical fault 
dimension (d = 2) HIHEl. 

In mean field theory, the local stress Tj on a given cell 
i is [13: 

n = J/Nj2( u 3- u i) + K L(vt-Ui) (4) 
3 

= Ju + K L vt - (K L + J)iii (5) 

where ut is the total offset of cell in the horizontal x di- 
rection, u — J^j U j/N is the average displacement, J/N 
is the mean-field elastic coupling strength between cells, 
and Kl ~ 1/VN is the loading stiffness 0] of the tec- 
tonic blocks that move far away from the fault with rela- 
tive velocity v. Initially the stresses Tj are randomly dis- 
tributed with T a i < Tj < t s j, where r s i is a fixed local 
static failure threshold stress and r a ^ is the fixed local ar- 
rest stress. The distributions of static failure stresses and 
arrest stresses represent the heterogeneity or geometrical 
disorder in the fault system. The differences between the 
failure and arrest stresses give the local distribution of 
stress drops during brittle failures; the earthquake dy- 
namics depend only on the stress drop distribution (no 
stochasticity). In addition, the scaling behavior of the 
system is not sensitive to the exact form of the distri- 
butions as long as they are bounded and T a ,j < T Sl j. 
We choose a compact distribution for r ai , such that 
p(r ai< ) = 3(W 2 - 4T a 2 ,)/(2VF 3 ) for -W/2 < t Qj1 < W/2 
and outside of these bounds. Also, we look at the low 
disorder limit where W <C T Sj j, where we choose T S) j = 1, 
so all cells will fail at this point. 

The fault is stuck while the stress on each cell is in- 
creased uniformly as dri/dt = K L v as a result of the 
external loading which is increased adiabatically (that 
is, we take the limit v — ► 0). When the stress on a 
cell reaches its failure threshold T S) j, the cell slips by the 
amount: 

Auj = (T s ,j - r a>i )/{K L + J). (6) 
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This stress drop is uniformly redistributed to all other 
cells (employed in the mean field approximation) by an 
amount: 

<5T J - = (c/A0(T., < -T Oli ),j^i (7) 

where c = J/(Kl + J) is the conservation parameter 
which gives the fraction of the stress drop of a cell that 
is retained in the system after it slips . The resulting 
stress increase on the other cells can cause some of them 
to slip as well, leading to an avalanche of cell slips, or an 
earthquake. 

DYNAMICAL WEAKENING AND 
STRENGTHENING 

The model includes dynamic weakening effects during 
the failure process |l7j, [l8|: after an initial slip in an 
earthquake, the strength of a failed cell is reduced to a 
dynamical value: 

Td,i = t s .i - e(r s ,i - T a .i), (8) 

with < e < 1 parameterizing the relative impor- 
tance of dynamical weakening effects in the system. This 
weakening represents the transition from static friction 
to dynamic friction during the rupture and the strength 
of a failed cell remains at the dynamic value throughout 
the remainder of the earthquake. In the time intervals be- 
tween earthquakes all failure thresholds heal back to their 
static value r Si i. Fisher et al. || found that at exactly 
e = the model produces a power law distribution of 
earthquake moments M following equation , cutoff by 
the finite fault size, with an analytical exponent (3 = 1/2 
(FIG. EJ- This corresponds to a b- value of 0.75, close to 
that associated with observed earthquakes on strike-slip 
faults 0. The power law scaling of the frequency-size 
statistics and other variables Q indicates that the model 
with e = operates at a critical point. In contrast, for a 
finite weakening e > the model produces the character- 
istic earthquake distribution, with power law statistics 
for the small events up to a cutoff moment that scales 
like: 

M c "* o// ~ 1/e 2 , (9) 

and quasi-periodically recurring large characteristic 
events that scale with the fault size (Mo ~ (LW) 3 / 2 ). 

The model can be expanded further to include dy- 
namic strengthening represented by e < 0. Ben-Zion 
and Sammis |'2-l| summarized multidisciplinary observa- 
tions which indicate that brittle failure of rock has an 
initial transient phase associated with strengthening, dis- 
tributed deformation, and creation of new structures. 
Detailed frictional studies also show an initial strength- 
ening phase associated with the creation of a new pop- 
ulation of asperity contacts 0,Hj|. In our model (FIG. 



Q we associate e < with regions off the main fault 
segments that are in an early deformation stage. To cap- 
ture basic aspects of brittle deformation on such regions 
in the three-dimensional volume around the main fault 
(FIG.(1J, we change the model as follows: when any cell i 
slips during an earthquake, and thereby reduces its stress 
by At; = Tf t i — r a ,i, the failure stress r/j of every cell 
j = 1, N is strengthenedby an amount |e| Ar</JV. Once 
the earthquake is complete, the failure stress of each cell 
is slowly lowered back to its original value. This repre- 
sents in a simple way the brittle deformation that occurs 
during an earthquake in the off- fault regions, which are 
first in a strengthening regime, compared to the main 
fault, and then have a weakening process. The events 
that are triggered as the failure stresses are lowered in 
the weakening period are referred to as aftershocks. The 
occurrence of aftershocks in this version of the model for 
off-fault regions is in agreement with the observation that 
a large fraction of observed aftershocks typically occur in 
off fault regions [2^| . For this version of the model with 
e < 0, both the primary earthquakes (i.e., main shocks) 
and the triggered aftershocks are distributed according 
to the Gutenberg- Richter distribution, up to a cutoff mo- 
ment scaling as 1/e 2 . Assuming that the increased fail- 
ure stress thresholds r/ ( j are slowly lowered with time as 
log(t) towards their earlier static values T Sl i, and that the 
stresses are distributed over a wide range of values, we 
show analytically in Appendix D that the temporal decay 
of aftershock rates at long times is proportional to 1 It, as 
in the modified Omori law of Eq. © with p= 1 0, M, m 
Remarkably, the long length scale behavior of this 
model can be shown to be the same as the behavior of the 
model given in Eq. JSJ with an added "antiferroelastic" 
term (— |e|Ju): 

n = Ju + K L vt - (K L + J)m - |e|Ju. (10) 

In Eq. (|10J) every time a cell fails, it slips by an amount 
Aitj that leads to stress loading of the other cells, lessened 
by |e| JAui/N compared to our original model (Eq. J5J). 
On the other hand, in the global strengthening model 
(described above) when a cell slips the failure stresses of 
all cells are strengthened by \e\JAm/N. On long length 
scales the global strengthening of the failure stress has 
equivalent effects on the earthquake statistics as the dis- 
sipation of the redistributed stress, up to corrections of 
order 0(1 /N), so the scaling behavior for large events of 
both models are the same. Moreover, Eq. HlOjl can be 
rewritten as: 

n = J[l - \e\][u - Ui ] + K L vt - [K L + J\e\]ui . (11) 

We can now absorb |e| by defining J' = J(l — |e|) and 
K' L = K L + J\e\. Rewriting Eq. ( |TT|> with the new 
definitions, and dropping the |e| contribution in [K' L — 
J\e\]vt since v — > 0, we find: 

n = J'u + K' L vt - {K' L + J')ui . (12) 
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Therefore we recover Eq. J5J with J — > J' and Kl — > 
A'^. This amounts to changing the stress conservation 
parameter c (from reference 19]). For Eq. 1)12(1 : 

c=J'/(K> L + J') = l-\e\ (13) 

where — > since we are concerned with the adiabatic 
limit. We also know (from reference |l9j) that the cut- 
off S c f for the Gutenberg-Richter distribution scales as 
S cf ~ 1/(1 - c) 2 . Thus, from Eq. lO we find that the 
cutoff for Eq. (HHJ) will scale as — l/|e| 2 . 

MAPPING TO SINGLE INTERFACE MAGNET 
MODEL 

The mean field version of the single interface magnet 
model with infinite range antiferromagnetic interactions 
is given bv [HE^: 

hi(t) = J[h - h t {t)} + H(t) ~kh + r)i(h) (14) 

where hi(t) is the position of the domain wall, H(t) is 
the external driving field, k is the coefficient of the an- 
tiferromagnetic term, and ru{h) is the pinning field. In 
the paper by Fisher et al. [8j it has been shown that the 
scaling behavior on long length scales resulting from Eq. 
((10JI . without the — |e| Ju term, is same as that of Eq. I|14(l 
without the antiferromagnetic term —kh. Furthermore, 
upon inspection we see the following correspondence be- 
tween the single interface magnet model (Eq. (|14|) ~). and 
the mean field earthquake model (Eq. (fTTH^ : 

- kh <S=S> -|e| Ju (15) 

In other words, the coefficient of the antiferromagnetic 
term k plays the same role in the magnet model (Eq. 
d}), as the coefficient of strengthening \e\J does in the 
earthquake model (Eq. (fLU|) '). 

PHASE DIAGRAM 

The regimes with various statistics produced by the 
model are summarized by the phase diagram given in 
FIG. EJ The range e > corresponds to "mature" lo- 
calized faults with a weakening rheology and character- 
istic earthquake statistics. The value e = corresponds 
to "immature" strongly inhomogeneous fault zones with 
Gutenberg-Richter statistics. Finally, the range e < 
corresponds to the fracture network away from the main 
fault, characterized by strengthening due to the creation 
of new structures and associated emerging aftershocks. 
It may be surprising that the discussed simple model can 
capture many of the essential general features of earth- 
quake statistics (or other systems with avalanches, such 



as driven magnetic domain walls). This can be under- 
stood through the renormalization group j^, |3(J , a pow- 
erful mathematical tool to coarse grain a system and 
extract its effective behavior on long space-time scales. 
Many microscopic details of a system are averaged out 
under coarse graining, and universal aspects of the be- 
havior on long scales depend only on a few basic prop- 
erties such as symmetries, dimensions, range of inter- 
actions, weakening/strengthening, etc. When a model 
correctly captures those basic features, the results pro- 
vide proper predictions for statistics, critical exponents, 
and universal scaling functions near the critical point. 
Consequently, many models that are in the same uni- 
versality class lead to the same statistics and exponents 
[1 03 Sill. The universal scaling functions around 
the critical point, discussed in the next section, provide 
additional information that can be used to distinguish 
between different models and universality classes. 

MOMENT RATE SHAPES 

In this section we focus on fault zones with Gutenberg- 
Richter power law statistics, modeled by systems at or 
close to the e = critical point. Recent analysis al- 
lowed researchers to obtain the moment rate dmo(t)/dt, 
which gives the slip on a fault per unit time during the 
propagation of earthquake rupture, for hundreds of large 
seismic events recorded on global networks |3lt |36j . The 
moment rates are derived from inversions of teleseismi- 
cally recorded seismograms on a global seismic network 
[sij ]. Motivated by works on statistical physics of mag- 
netic systems discussed in the previous chapter (see also 
[23ll30l p. we are interested in studying the event-averaged 
moment rate time profile (FIG. EJ) for earthquakes with 
given total moment Mq, denoted with (dmo(t\Mo)/dt) 
and the event-averaged moment rate time profile (FIG. 
0J) for earthquakes with given duration T, denoted with 
(dmo(t\T)/dt). Here mo(t\T) is the (cumulative) mo- 
ment at time t of the propagating earthquake of total 
duration T, and rrio(t\Mo) is the cumulative moment at 
time t of the earthquake of total moment Mo. Theoret- 
ical analysis of phase diagrams similar to that shown in 
FIG. H implies that near the critical point there should 
be, in addition to scaling exponents, also universal scaling 
functions (up to a rescaling of the ordinate and abscissa) 
[23j. In our model the two scalable functions of interest, 
(dmo(t\M )/dt} and (dmn(t\T) I dt ) , obey respectively the 
following scaling relations |32l l33| : 

(dm (t\M )/dt)/A4 /2 ~ /(</M 1/2 ) (16) 

and 

{dm (t\T)/dt) ~ g(t/T) (17) 
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We determined these scaling functions from corre- 
sponding results for magnets, using the fact that our 
mean field version of the Ben-Zion and Rice model of 
Eq. (JjJJ 8, 19] is in the same universality class (i.e. has 
the same universal behavior on long length scales) as the 
above mentioned model for domain wall motion in mag- 
nets. 



try may result from a rupture process that begins with 
a failure of a large asperity, from finite fault size effects 
if the rupture process slows down once the rupture has 
traversed the fault in one direction, or from contributions 
of early aftershocks. 

DISCUSSION 



EXPONENTS AND DATA COLLAPSES 

We compare the observation results with our model 
and find remarkable agreement in most cases. The 
frequency-moment distribution, D{Mq) ~ Mq 1 ^^ of the 
observed data [sH has (inset of FIG.I2J three decades of 
scaling and an exponent of /3 = l/2±0.05, in close agree- 
ment with the model near e = 0. The deviation from a 
power law distribution at the low moment range is as- 
sociated with the reduced resolution of the observational 
network for small events. In mean field theory the univer- 
sal scaling function f(x) in Eq. (|16|) is of the exact form 
[33|: f m f(x\ = Axe~ Bx I 2 with non-universal constants 
A = B = 1. In FIG.Olwe present a collapse of the obser- 
vational data of (dmo(t\AIo)/dt) for four different values 
of Mq to obtain the corresponding function f exp {x) for 

1/2 

observations with x = t/M . The observational curves 
not only collapse, and are therefore universal, the mean 
field exponent 1/2 in the scaling variable x is in excel- 
lent agreement with observations. We fit the functional 
form f m f(x) with A = 4 and B = 4.9 to the collapse of 
the observed data; f eX p(x) deviates from the f m f(x) for 
small values of the ordinate. 

In mean field theory, the function g(x) of Eq. I|17|l 
is of the symmetric form: g m f{x) — Ax{\ — x), where 
A is a non-universal constant. In FIG. 0] we collapse 
observational data for (dmo(t\T)/dt) for three values of 
T to obtain the function g exp (x) with x = t/T. Again we 
find that the curve collapses well, even though only small 
data sets were available, and the exponent of 1 obtained 
from the collapse is in excellent agreement with mean 
field theory. 

We also plot the mean field scaling function g m f(x) 
with A = 80. The results show that while the scaling ex- 
ponents agree, there are notable differences between the 
observational function g exp {x) and the mean field func- 
tion g m f(x), especially for small values of the ordinate. 
We checked that finite size effects do not play a role in 
9m fix) (or in f m f{x) for that matter). Also, we find 
that the mean skewness coefficient for the g exp {x) curves 
is 0.878 and the mean standard error of skewness is 0.705: 
since twice the standard error is greater than the abso- 
lute value of the skewness, the asymmetry is not sta- 
tistically significant. Therefore more observational work 
with a larger data set is required to verify the moment 
rate shape asymmetry and clarify its origin. An asymme- 



The employed earthquake model |l7j, [18| was shown 
in the past to have a critical point at e = and addi- 
tional dynamic regimes for e>o|l 0] compatible with 
observed frequency-size statistics of earthquakes on indi- 
vidual fault zones 0, Q] . We have generalized the the- 
oretical analysis to include a strengthening regime e < 
with aftershocks, and derived universal scaling functions 
around the critical point e — 0. The results provide new 
tools for data analysis that may be used to obtain an 
improved understanding of earthquake dynamics. The 
analysis indicates that near e = the model is in the 
same universality class as a recent model for domain wall 
motion in magnets and we can match the phase diagram 
regions e < and e > to those of the corresponding- 
magnet model I n other words, the two sys- 
tems are marked by the exact same universal scaling ex- 
ponents, universal scaling functions, and similar phase di- 
agrams. The model predictions for frequency-size statis- 
tics and moment-rates of earthquakes near e = are 
overall in close agreement with observational data of rel- 
atively large earthquakes recorded on the global seismic 
network |3l| . However, the observed mean moment-rate 
of earthquakes with a given duration T apparently in- 
creases with time more rapidly than it drops off, con- 
trary to the corresponding symmetric model function. 
The potential asymmetry in the observed data may re- 
sult from a rupture process that begins with a failure of 
a large asperity. This is compatible with observations 
that hypocenter locations tend to be located close to an 
area on a fault that produces large moment release, (e.g. 
[IH). As we explain below, an analogy to magnetic sys- 
tems suggests that an asymmetry could potentially also 
stem from momentary initial threshold strengthening as- 
sociated with the creation of a new population of asper- 
ity contacts upon local failure and followup aftershocks. 
Further study will be necessary to clarify this issue. 

Theoretical analysis of such a potential asymmet- 
ric rupture process requires corrections to the mean 
field earthquake model results. Recently it has been 
shown that the corresponding magnetic domain-wall 
model [13, H3 predicts well the critical scaling exponents 
for Barkhausen noise experiments in magnets. Signifi- 
cantly, the experimental scaling function for magnetiza- 
tion avalanches or Barkhausen "pulses" HSHil; that is 
the analogue of the moment-rate time profile for fixed 
earthquake duration of Eq. I|17|) , shows the same type of 
asymmetry that is apparently observed for earthquakes 
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(FIG. 0}. It has been suggested that the asymmetry 
in the function is due to eddy currents in the magnet 
[35| . Eddy currents have a similar effect in magnets as 
transient threshold strengthening would have on earth- 
quakes. In this paper we have shown how long-term 
threshold strengthening (on time scales longer than in- 
dividual earthquakes) leads to aftershocks, which effec- 
tively represent an asymmetry on time scales longer than 
individual earthquakes. This raises the possibility that 
the origin of this asymmetry may be similar in both 
magnets and earthquakes: in both cases it may be due 
to a transient force (due to eddy currents or threshold 
strengthening respectively) that counteracts the propa- 
gation of an event and thus leads to asymmetric event 
profiles that taper off more slowly than they began. Our 
study shows that there are important theoretical and ob- 
servational connections between processes in earthquake 
and magnetic systems. 

A related study of earthquake moment rate shapes was 
done by Houston [3(| who used a data set similar to the 
data of Bilek [3l| used here. In order to compare the av- 
erage moment rate profiles of these earthquakes, the pro- 
files were rescaled by two methods: moment scaling and 
duration scaling. Moment scaling rescales the time axis 
of the profiles using assumptions about cracklike scaling 
of the events (different from our mean- field result) , and 
then rescales the height of the profiles so that the area 
(or moment) under all profiles are the same. The dura- 
tion scaling of [jjfj rescales the time axis such that all 
duration-scaled profiles end at the same reference dura- 
tion, and then rescales the moment rates so that all scaled 
time profiles have the same area underneath. 

Fig. El shows average moment rate profiles that 
were obtained by Houston from moment scaling (top) 
and duration scaling (bottom) for data from sev- 
eral subduction zones. The top plots correspond to 
(dmo(t\Mo) I dt) collapsed, and the bottom plots corre- 
spond to (dmo(t\T)/dt) collapsed. We see from the top 
part in Fig. 0that overall the moment rate shapes seem 
to agree rather well with our predicted mean field shapes 
of Fig. |3J Likewise the bottom plots of Fig. appear to 
agree quite well the mean field curve of Fig. 0] except for 
an additional slight asymmetry, with positive skewness, 
similar to our result from observations. 

In [3^| the skewness of the duration scaled moment rate 
profiles, i.e. the skewness of the (dmo(t\T) / dt) collapses, 
is calculated and found to range from 0.12 to 0.36 depend- 
ing on the depth of the earthquakes (bigger skewness for 
greater depth) and the method used to extract moment 
rates from seismograms. Since the statistical error of the 
skewness is not given in [3(|, we cannot determine if the 
above skewness values are statistically significant. It is 
interesting, however, that just like in the case of mag- 
nets, and in our analysis of the observational data, the 
skewness is slightly positive, indicating that on average 
earthquakes tend to grow faster than they die down. The 



skewness values quoted in [36J are smaller than the 0.878 
value obtained here, though one of them (the 0.36 value 
for earthquakes at greater depth) is within errorbars. 

While overall the shapes of the moment rate profiles 
look rather similar for both Houston's analysis and our 
mean field theory results, there are significant differences 
in the analysis that can be summarized as follows: 

(1) The exponents used in [36| to create the moment- 
scaled average profile (dmo(t\Mo) / dt) differ from the ex- 
ponents we used. In both cases the x axis is divided 
by Mq* and the y axis is divided by M^~ at in order to 
ensure normalization of the curves to a fixed reference 
moment. In pjjj the exponent value is at — 1/3, while in 
this paper we use the mean field value at — 1/2 (see Fig. 
EJ. The value a t = 1/3 in [H is associated with crack- 
like scaling; for cracks the moment Mg scales with the 
rupture area A as Mq ~ A 3 / 2 4] . Since the rupture area 
scales with the diameter L of the crack as A ~ I? and 

the duration T of the earthquake scales with its diameter 

1/3 

as T ~ L one arrives at the proposed scaling T ~ M . 
In [3(| it is mentioned, however, that the data would also 
be compatible with rescaling the time (or x) axis with an 
exponent of at — 0.41, i.e. that T ~ Mq' 41 for the data. 
This scaling is much closer to the mean field prediction 

1 /2 

T ~ M ' , i.e. at = 1/2, that is used in this work. 

Our exponent at — 1/2 is derived from the mean field 
prediction || that the duration T of earthquakes in the 
critical power law region of Fig. 13 scale as T ~ M Q 1/2 . In 
8] it is shown analytically that the moment Mo of earth- 
quakes in the critical regime scales with the rupture area 
as Mo ~ A. Since the rupture area scales, as before, with 
the earthquake diameter as A ~ L 2 , and the earthquake 
duration scales as T ~ L, one obtains for the mean field 

1 /2 

prediction T ~ M Q for critical earthquakes. Note that 
for earthquakes which are so large that their horizontal 
diameter L is much larger than the vertical width of the 
fault, one expects the moment to scale as Mo ~ L. Since 
again the duration scales as T ~ L, one obtains T ~ Mo, 
i.e. at = 1. The data used here and in 36] apparently 
does not fall within this scaling regime. 

(2) Our theoretical analysis leads to predictions of both 
scaling exponents and entire scaling functions, which can 
be compared to data. In contrast, in [3j| the shape is 
obtained only empirically by averaging rescaled data, and 
it is dependent on the chosen scaling exponents. Our 
procedure of comparing the theoretical results to data 
is designed to provide simultaneously, through a scaling 
collapse that minimizes deviations of individual averaged 
curves from the collapsed one, the scaling exponent at 
and the entire scaling function. On the other hand, in 
[3^| the value of at is obtained from a separate analysis of 
duration versus moment of observed data, and the shape 
is then obtained empirically, by averaging rescaled data. 

(3) For our simulated moment rates, the exponents 
used in |36J would not lead to a good collapse, i.e. the 
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deviation between various rescaled curves would be larger 
for at = 1/3 compared to a t = 1/2. 

Unfortunately, since the available observational data 
have large error bars, we cannot determine at present un- 
equivocally, which set of exponents works better. More 
moment rate data, especially for many small earthquakes 
would allow us to reduce statistical errorbars and obtain 
a more precise comparison of the scaling functions in the- 
ory and experiment. The smaller earthquakes are also 
more likely to fall into the critical scaling regime of our 
mean field prediction. However, moment rate data are 
more difficult to obtain for small earthquakes, as they re- 
quire greater spacial resolution. We hope that this study 
will motivate more observational work and analysis to 
answer these questions. 



DETERMINATION OF OMORI LAW OF 
AFTERSHOCK DECAY 

We derive the Omori law for aftershock decay from 
our mean field earthquake model in the case of dynamic 
strengthening. After a primary earthquake, the static 
failure stress, r s , is increased by an amount e(s) /N, where 
(s) is the mean event size. Given in terms of moment, 
the mean avalanche size in mean field theory is simply 
given by (l9| : 



1A 



s/s 3/2 ds ~ 1/c 



(18) 



where 1/s 3 / 2 ~ D(s), is the scaling behavior of the mean 
field event size distribution. Now since (s) ~ 1/e, the 
failure stress, r s , is then shifted by an amount 1/iV. 

In order to produce aftershocks, we assume that the 
new larger failure stress, which we call t2, is lowered 
slowly, logarithmically in time until it returns to its orig- 
1/N. More precisely, using t = 



inal value, r. 
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for the time at the end of the previous earthquake, t 
decays as follows: 



/ 



T f (t)=r f (l-log(l+t)) 



(19) 



We define rj(t) as the amount by which t - has been low- 
ered at time t: 



n{t)=T° f -T f (t)=T° f log(l + t) 



(20) 



Eq. Ijl9|l is only accurate to first order since lowering the 
threshold triggers on average ~ r](t) aftershocks (assum- 
ing nonsingular stress distributions). Since aftershocks 
occur quickly compared to the weakening time scale, the 
aftershocks that are triggered due to lowering the thresh- 
old will cause T/(i) to shift up by: T](t) x e(s)/N ~ 
r)(t)/N. Thus the expression for Tf(t) becomes: 

r f (t) = 7-0(1 - log{\ + t)) + n(t)/N (21) 



and Eq. J3UJ gives: 

n(t)=r° f log(l + t)- V (t)/N 

r,{t) 



TjlogQ- + t) 



l + l/N 



(22) 
(23) 



If we assume that the number of aftershocks triggered by 
time t is A{t) ~ we obtain the modified Omori law 
(Eq.©) for aftershock decay 



dA(t) 
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(24) 



dt (l + t)(l + l/N) 
For large t and large N we have the original Omori law 
dA(t) r° 



dt 



t 



(25) 
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FIG. 1: A planar representation of a 3-D segmented fault 
by a 2-D heterogeneous fault embedded in a 3-D solid 
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FIG. 2: Phase diagram of the model (see text and |19| for 
details) . 
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FIG. 3: A collapse of averaged earthquake pulse shapes, 
(dmo(t\Mo) I dt) , with the size of the moment Mo in Newton 
meters within 10% of each size given in the legend respec- 
tively. In order to obtain each collapsed moment rate shape, 
five to ten earthquakes were averaged for each value of Mo- 
The collapse was obtained using the mean field scaling re- 
lation H: (drn (t\M )/dt}/Mf/ 2 ~ f(t/M 1/2 ) (see text Eq. 

In our mean field theory the universal scaling func- 
tion is f m f(x) = Axe~ Bx2/2 where x = t/M Q 1/2 . We plot this 
functional form (bold curve) with A = 4 and B — 4.9. Inset: 
The raw data and the averaged data (before collapsed). 
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FIG. 4: A collapse of averaged earthquake pulse shapes, 
(dmo(t\Mo) /dt) with a duration of T (seconds) within 10% 
(given in legend), is shown. The collapse was obtained using 
the mean field scaling relation H3: (dmo(t\T)/dt) ~ g(t/T) 
In order to obtain each collapsed pulse shape, two to 
ten earthquakes were averaged for each value of T. In our 
mean field theory the universal scaling function is g m f{x) = 
Ax(l — x) with x = t/T. We plot this functional form (bold 
curve) with A = 80. Note the apparent asymmetry to the left 
in the observed data while the theoretical curve is symmetric 
around its maximum. Inset: The raw data and the averaged 
data (before collapsed). 
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FIG. 5: Reseated moment rate versus time profiles of seis- 
mic events determined with moment scaling (top) and dura- 
tion scaling (bottom) obtained from |3r^| (see text). The data 
shown are obtained from several subduction zones in the indi- 
cated geographical locations and the numbers of earthquakes 
used for each region are given in parentheses. 



